#include "Twogrid.h"
#include <fstream>

//1D
double f1(double *x);
double g1(double *x);
//2D
double f2(double *x);
double g2(double *x);

int main(int argc, char *argv[])
{
    TwoGrid2D Solver;
    std::vector<std::vector<double>> Region(2);
    Region[0].resize(2);
    Region[0][0] = 0.0;
    Region[0][1] = 0.0;
    Region[1].resize(2);
    Region[1][0] = 1.0;
    Region[1][1] = 1.0;
    int n = 49;
    std::valarray<double> init(0.0, n);
    Solver.Initialize(init, f2, g2, Region);
    std::string RS = "Weighted_Jacobi";
    double w = 1.0;
    Solver.Set_relaxation_scheme(RS, w);
    double a = 1e-6;
    int max = 20;
    Solver.Set_termination_conditions(a, max);
    std::valarray<double> solution(n);
    solution = Solver.Solve();

    int m = (int)sqrt(n);
    double h = 1.0/(m+1);
    std::ofstream output("project.m");
    output << "x = zeros(1," << n << ");" << std::endl;
    output << "y = zeros(1," << n << ");" << std::endl;
    output << "u = zeros(1," << n << ");" << std::endl;
    for(int i = 0; i < m; i++)
    {
	for(int j = 0; j < m; j++)
	{
	    output << "x(" << j+i*m+1 << ") = " << (j+1)*h << ";" << std::endl;
	    output << "y(" << j+i*m+1 << ") = " << (i+1)*h << ";" << std::endl;
	    output << "u(" << j+i*m+1 << ") = " << solution[j+i*m] << ";" << std::endl;
	}
    }
    output << "T = delaunay(x,y);" << std::endl;
    output << "trisurf(T,x,y,u);";
    output.close();

    return 0;
};

double f1(double *x)
{
    //return 2;
    return sin(x[0]);
};

double g1(double *x)
{
    //return x[0]-x[0]*x[0];
    return sin(x[0]);
};

double f2(double *x)
{
    return sin(x[0])*sin(x[1]);
};

double g2(double *x)
{
    return 2.0*sin(x[0])*sin(x[1]);
};



